home *** CD-ROM | disk | FTP | other *** search
/ HPAVC / HPAVC CD-ROM.iso / MOR55SRC.ZIP / MORIA / SOURCE / RND.C < prev    next >
C/C++ Source or Header  |  1992-12-07  |  3KB  |  117 lines
  1. /* source/rnd.c: random number generator
  2.  
  3.    Copyright (c) 1989-92 James E. Wilson
  4.  
  5.    This software may be copied and distributed for educational, research, and
  6.    not for profit purposes provided that this copyright and statement are
  7.    included in all such copies. */
  8.  
  9. #include "config.h"
  10. #include "constant.h"
  11. #include "types.h"
  12.  
  13. /* Define this to compile as a standalone test */
  14. /* #define TEST_RNG */
  15.  
  16. /* This alg uses a prime modulus multiplicative congruential generator
  17.    (PMMLCG), also known as a Lehmer Grammer, which satisfies the following
  18.    properties
  19.  
  20.    (i)     modulus: m - a large prime integer
  21.    (ii)     multiplier: a - an integer in the range 2, 3, ..., m - 1
  22.    (iii) z[n+1] = f(z[n]), for n = 1, 2, ...
  23.    (iv)     f(z) = az mod m
  24.    (v)     u[n] = z[n] / m, for n = 1, 2, ...
  25.  
  26.    The sequence of z's must be initialized by choosing an initial seed
  27.    z[1] from the range 1, 2, ..., m - 1.  The sequence of z's is a pseudo-
  28.    random sequence drawn without replacement from the set 1, 2, ..., m - 1.
  29.    The u's form a psuedo-random sequence of real numbers between (but not
  30.    including) 0 and 1.
  31.  
  32.    Schrage's method is used to compute the sequence of z's.
  33.    Let m = aq + r, where q = m div a, and r = m mod a.
  34.    Then f(z) = az mod m = az - m * (az div m) =
  35.          = gamma(z) + m * delta(z)
  36.    Where gamma(z) = a(z mod q) - r(z div q)
  37.    and     delta(z) = (z div q) - (az div m)
  38.  
  39.    If r < q, then for all z in 1, 2, ..., m - 1:
  40.    (1) delta(z) is either 0 or 1
  41.    (2) both a(z mod q) and r(z div q) are in 0, 1, ..., m - 1
  42.    (3) absolute value of gamma(z) <= m - 1
  43.    (4) delta(z) = 1 iff gamma(z) < 0
  44.  
  45.    Hence each value of z can be computed exactly without overflow as long
  46.    as m can be represented as an integer.
  47.  */
  48.  
  49. /* a good random number generator, correct on any machine with 32 bit
  50.    integers, this algorithm is from:
  51.  
  52. Stephen K. Park and Keith W. Miller, "Random Number Generators:
  53.     Good ones are hard to find", Communications of the ACM, October 1988,
  54.     vol 31, number 10, pp. 1192-1201.
  55.  
  56.    If this algorithm is implemented correctly, then if z[1] = 1, then
  57.    z[10001] will equal 1043618065
  58.  
  59.    Has a full period of 2^31 - 1.
  60.    Returns integers in the range 1 to 2^31-1.
  61.  */
  62.  
  63. #define RNG_M 2147483647L  /* m = 2^31 - 1 */
  64. #define RNG_A 16807L
  65. #define RNG_Q 127773L       /* m div a */
  66. #define RNG_R 2836L       /* m mod a */
  67.  
  68. /* 32 bit seed */
  69. static int32u rnd_seed;
  70.  
  71. int32u get_rnd_seed ()
  72. {
  73.   return rnd_seed;
  74. }
  75.  
  76. void set_rnd_seed (seedval)
  77. int32u seedval;
  78. {
  79.   /* set seed to value between 1 and m-1 */
  80.  
  81.   rnd_seed = (seedval % (RNG_M - 1)) + 1;
  82. }
  83.  
  84. /* returns a pseudo-random number from set 1, 2, ..., RNG_M - 1 */
  85. int32 rnd ()
  86. {
  87.   register long low, high, test;
  88.  
  89.   high = rnd_seed / RNG_Q;
  90.   low = rnd_seed % RNG_Q;
  91.   test = RNG_A * low - RNG_R * high;
  92.   if (test > 0)
  93.     rnd_seed = test;
  94.   else
  95.     rnd_seed = test + RNG_M;
  96.   return rnd_seed;
  97. }
  98.  
  99. #ifdef TEST_RNG
  100.  
  101. main ()
  102. {
  103.   long i, random;
  104.  
  105.   set_rnd_seed (0L);
  106.  
  107.   for (i = 1; i < 10000; i++)
  108.     (void) rnd ();
  109.  
  110.   random = rnd ();
  111.   printf ("z[10001] = %ld, should be 1043618065\n", random);
  112.   if (random == 1043618065L)
  113.     printf ("success!!!\n");
  114. }
  115.  
  116. #endif
  117.